# the factorial function
def fak n
  return 1 if n==1 || n==0
  n*fak(n-1)
end

# returns the product of the numbers in the range
def pi(from,to)
  sum=1
  (from+1..to).each{|j|
    sum*=j
  }
  sum
end

#returns the binomial coefficient.
#the number of ways that B objects can be chosen from A objects.
def binom(a, b)
  fak(a)/(fak(b)*(fak(a-b)))
end


#returns binom(cards-aces,draws-drawAces) / binom(cards,draws).
#it's not calculated in the straightforward way because then we would get too big numbers in intermediate results.
#instead we have moved a bit around with the formula.
def f(aces,drawAces,cards,draws)
  a=1.0/(pi(cards-aces,cards))
  b=pi(draws-drawAces,draws)
  c=pi(cards-aces-draws+drawAces,cards-draws)
  a*b*c
end

#returns the probability of drawing exactly drawAces.
#which is ( binom(aces,drawAces)*binom(cards-aces,draws-drawAces) )/ binom(cards,draws).
def propExact(aces,drawAces,cards,draws)
    binom(aces,drawAces)*f(aces,drawAces,cards,draws)
end

#returns the probablility of drawing at least drawAces aces
def prop(aces,drawAces,cards,draws)
  sum=0
  (drawAces..aces).each{|i|
    sum+=propExact(aces,i,cards,draws)
  }
  sum
end

#example calculation:
#puts the probability of drawing at least 3 aces when you draw 65000 cards.
#The deck of cards has 650000 cards and 5 aces.
puts prop(17,14,650000,65000)


